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ABSTRACT 

We present results of numerical simulations carried out with a 2D radiation 
hydrodynamics code in order to study the impact of massive stars on their sur- 
rounding interstellar medium. This first paper deals with the evolution of the 
circumstellar gas around an isolated 60 star. The interaction of the pho- 
toionized H ll region with the stellar wind bubble forms a variety of interesting 
structures like shells, clouds, fingers, and spokes. These results demonstrate that 
complex structures found in H ll regions are not necessarily relics from the time 
before the gas became ionized but may result from dynamical processes during 
the course of the H ll region evolution. We have also analyzed the transfer and 
deposit of the stellar wind and radiation energy into the circumstellar medium 
until the star explodes as a supernova. Although the total mechanical wind en- 
ergy supplied by the star is negligible compared to the accumulated energy of the 
Lyman continuum photons, the kinetic energy imparted to the circumstellar gas 
over the star's lifetime is 4 times higher than for a comparable windless simula- 
tion. Furthermore, the thermal energy of warm photoionized gas is lower by some 
55 %. Our results document the necessity to consider both ionizing radiation and 
stellar winds for an appropriate description of the interaction of OB stars with 
their circumstellar environment. 

Subject headings: galaxies: evolution — H ll regions — hydrodynamics — instabilities — 
ISM: bubbles — ISM: structure 
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1. Introduction 

Massive stars play an important role in the evolutionary history of galaxies. They are 
the primary source of metals and dominate the turbulent energy input into the interstellar 
medium (ISM) by stellar winds, radiation, and supernova explosions. The radiation field of 
a massive star first dissociates the ambient molecular gas and forms a so-called photodis- 
sociation region of neutral hydrogen. Subsequently, the Lyman continuum photons of the 
star ionize the H I gas and produce an H ll region that expands into the neutral ambient 
medium. A fast stellar wind creates shocks that form a so-called stellar wind bubble (SWB) 
filled with hot plasma, which expands into the H ll region. Finally, the star explodes as a 
supernova of typell (SN ll), creating a supernova remnant (SNR) that sweeps up the ambi- 
ent medium. The SWBs and SNRs of neighboring stars can overlap and form a superbubble 
with a diameter of order 1 kpc. 

This paper pursues two major goals. First, we want to examine the combined influence 
of wind and ionizing radiation on the dynamical evolution of circumstellar matter around 
massive stars, i.e. we are interested in the interaction processes between the photoionized 
H II region and the SWB that evolves into the ionized gas. The second goal is to improve 
our knowledge of the energy transfer efficiency between massive stars and the ISM: How and 
to what fraction is the energy of stellar radiation and the stellar wind converted into kinetic, 
thermal, and ionization energy of the ISM? How does the formation of the SWB infiuence 
the transfer of stellar radiation? To what extent does all this depend on the evolutionary 
state of the star? 

To investigate these effects we perform numerical 2D radiation hydrodynamic simula- 
tions of the interaction between an isolated massive star and its surrounding ISM via stellar 
hydrogen-ionizing photons and a stellar wind. We calculate the hydrodynamical evolution 
of the circumstellar gas coupled with radiation transfer, time-dependent ionization of hydro- 
gen, and a realistic description of cooling. The stellar mass-loss rate, the terminal velocity 
of the wind, the effective temperature, and the luminosity of the star are specified as time- 
dependent boundary conditions. We examine the evolution of the circumstellar material 
starting from the main sequence (MS) phase of the star until it explodes as a supernova. We 
do not consider the SNR formation. In this paper we present results of these calculations 
for a star with an initial mass of 60 Mq, extending the work of Garcia-Segura, Mac Low, & 
Langer (1996a, hereafter GMLl) to a more precise analysis of the energetic aspect during 
the whole evolution of the star. 

The remainder of this paper is structured as follows: In section 2 we briefly review 
theoretical and observational studies of H ll regions and SWBs around massive stars for later 
comparison. Our numerical procedure as well as the initial conditions and time-dependent 
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boundary conditions are described in section 3. In section 4 we present the numerical results 
of our calculations and discuss them in the context of previous analytical approaches, other 
numerical investigations, and observations. Finally, in section 5 we summarize our main 
results and draw some conclusions from the work presented in this paper. 



2. Theoretical concepts and observational constraints 

2.1. H II regions 

H II regions are observed in various sizes and shapes. Depending on their geometrical 
extent, they are labeled as ultracompact, compact, extended, or giant. The ultracompact 
(linear size < 0.1 pc) and partly also the compact (0.1-0.3 pc) H ll regions are still deeply 
embedded in their birthplaces, mature molecular clouds, and due to dust obscuration gen- 
erally are only observable in the IR and radio wavelengths. On the other hand, extended 
H II regions with sizes up to several parsecs and giant H ll regions which are composed of 
individual H ll regions and have sizes of a few 100 pc are often easily discernible in the optical 
by their bright Balmer and forbidden metal lines. Large holes and shells make it obvious 
that the giant H II regions are not only powered by stellar radiation but also by stellar winds 
and supernova explosions in the underlying OB star cluster (see e.g. Yang et al. 1996). 

Within the framework of the simplest theoretical approach an O star suddenly "turns 
on" in a constant density, motionless medium and begins to ionize its surroundings. The 
reaction of the medium to the stellar photons is well-known and has been described in detail 
in standard textbooks (see e.g. Spitzer 1978). For comparison purposes with our numerical 
results we are particularly interested in the time evolution of the location ri of the ionization 
front: 

n = Qcs^nro^V + rJ/^^ (1) 
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where Cg^n is the isothermal sound speed of the ionized gas, Tq the initial Stromgren radius, 
LhyC the photon luminosity in the Lyman continuum, no the hydrogen number density in the 
neutral ambient medium, /^^^^ (T) the coefficient for recombinations of hydrogen into levels 2 
and higher, and r the age of the star. According to Lasker (1967) the kinetic energy of the 
expanding swept-up shell is 
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where mn is the mass of a hydrogen atom. Because this swept-up shell can be expected to 
remain cool and neutral due to strong radiative cooling, we can also estimate the ionization 
energy stored in the H ll region: 



Here we have used equation (2) for the time-dcpcndcnt radius of the ionization front and Xo 
is the ionization potential of hydrogen in the ground state. 

The same ansatz can then be used to estimate the thermal energy of warm gas in the 
H II region: 



where we will insert Tn = 8000 K as an appropriate approximation to the temperature in 
the H II region, /cb is the Boltzmann constant. 



First considerations about the inner structure of SWBs go back to Pikel'ner (1968), 
Avedisova (1972), Dyson & de Vries (1972), and Dyson (1973), but the work of Castor, 
McCray, & Weaver (1975) and Weaver et al. (1977) set a milestone in the approach of 
understanding SWBs. They presented a fairly complete picture of the structure and evolution 
of SWBs together with a set of equations that describes the evolution under the simplifying 
assumptions of a point source of a constant and spherically symmetric strong wind that 
interacts with a homogeneous ambient ISM. The global structure that arises from such a 
wind-ISM interaction is depicted in Figure 1 together with the expanding H II region into 
which the SWB evolves. It consists of a free-flowing supersonic wind that is heated to about 
10^ — 10^ K when it passes the inner reverse shock at Tgi- The pressure of this hot rarefled gas 
that normally fills most of the volume of the SWB is typically much higher than the pressure 
in the photoionized ambient medium. As a consequence, the hot gas bubble expands into 
the H II region producing a forward shock at that sweeps up the gas from the H ll region 
in a shell which is separated from the hot bubble interior by a contact discontinuity at Tc. 

Three phases in the evolution of such a bubble can be distinguished: The first is the 
adiabatic phase that lasts until the shock speed Vs2 drops below 200kms~^ (Falle 1975), 
typically a few 10 to a few 1000 yr, depending on the mechanical luminosity of the wind and 
on the ambient density. The bubble is expanding so fast that radiative cooling does not play 
a significant role for the dynamical behavior. The second stage of evolution is characterized 







2.2. Theory of stellar wind bubbles 
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by strong cooling in the shell of swept- up material (between Tc and rs2), allowing it to be 
compressed into a geometrically thin (rdt) = rs2{t), where t denotes time), dense shell, 
whereas for the hot bubble (between r^i and r^) cooling is still negligible. This phase lasts 
much longer than the first, so-called adiabatic phase. Thermal conduction from the hot 
bubble interior to the collapsed shell may become important and may modify the structure 
of the bubble in this stage. Weaver et al. (1977) assume that an equilibrium between the 
conductive energy flux and the mechanical energy flux due to evaporation of shell mass in 
the reverse direction is established. Under the assumptions that in the hot bubble between 
Tsi and Tc = rs2 the pressure is everywhere the same and that the thermal energy contained 
within this region is much higher than the kinetic energy. Weaver et al. (1977) derived the 
following equations for the temporal dependence of the radius rs2{t) and the pressure in the 
hot bubble Ph{t) in this stage: 

rAt) = = (l^j'^ii'Vo'^t^'S (6) 

Po = no mu is the density of the ambient medium and = ^M^v^ the stellar wind lumi- 
nosity, with Mw and Vw being the mass-loss rate and the terminal velocity of the star. 

It is also possible to obtain an analytical result for the energy that is deposited in the 
ISM by SWBs. We neglect the first (fully adiabatic) stage of bubble evolution because it is 
very short. Cooling in the hot bubble strongly depends on the efficiency of heat conduction, 
or, how efficiently heat conduction is suppressed by magnetic fields. We will assume for 
the analytical theory that there is no heat conduction and thus cooling in the hot bubble is 
unimportant during the lifetime of the star. The hot bubble acts like a piston on the ambient 
medium and analogous to the case of an H ll region it can be shown that the PdV work is 
equally distributed on kinetic energy of shell motion and thermal energy in the shell (which 
gets immediately lost due to cooling). If we furthermore assume that the kinetic energy 
of the stellar wind is completely transformed into thermal energy at the reverse shock (the 
strong-jump conditions give 15/16, but dumpiness of the stellar wind, which we completely 
neglect, might reduce this value), we get for the kinetic energy: 



1 



Vs2{t) 



£^k=2 / P^^Vs2. (8) 







Vs2 is the volume inside rs2- Using equation (6) for rs2 and equation (7) for the pressure in 
the hot bubble yields 

E^^^L^T. (9) 



-6- 



The thermal energy of the hot gas is thus 

P,>dV,2 = —L^T . (10) 

Besides heat conduction and coohng in the hot bubble we have also neglected the fact that 
part of the thermal energy might be used for coUisional ionization of gas. Thus, the results 
for £'k and Et are upper limits. 

Though the analytical and semi-analytical solutions for the evolution of SWBs have been 
improved in recent years (Koo & McKee 1992a,b; Garcfa-Segura & Mac Low 1995; Pittard, 
Dyson, & Hartquist 2001a; Pittard, Hartquist, & Dyson 2001b) a variety of physical effects 
remains to be included in order to achieve a better agreement of models and observations. 
For example, the discrepancy between models and observations with regard to the evolution 
of the hot phase in bubbles has recently been reviewed by Mac Low (2000). See also Chu 
(2000). It has become evident that the stellar parameters such as the mass-loss rate, the 
terminal velocity, the effective temperature and the luminosity of the star which drive the 
evolution of the circumstellar matter vary strongly over time. Because most previous studies 
dealt with either the evolution of H II regions or of SWBs, little is known about the interaction 
of these two structures. 



2.3. Observations of stellar wind bubbles 

Hot gas (10^ — 10^ K) is expected to be a crucial indicator for the existence of SWBs. 

However, there have been only two successful X-ray observations of SWBs so far: NGC 6888 
(Bochkarev 1988; Wriggc, Wendker, & Wisotzki 1994) and S308 (Wrigge 1999). Both are 
actually Wolf-Rayet (WR) bubbles (i.e. the wind-driving star has already entered its WR 
stage). The fact that no MS bubble has yet been observed in X-rays may be due to their 
physical properties; they are expected to be large, diffuse, and consequently dim (compared 
with WR bubbles). There may also be a problem with the theoretical models. More sensitive 
X-ray telescopes of the next generation should be able to illuminate this problem. 

The shells of MS bubbles are also difficult to observe in the optical because they are 
large and dim (McKee, Van Buren, & Lazareff 1984). Nevertheless, shells around MS stars 
as well as shells around evolved massive stars, which — based on their radius, expansion 
velocities, and shell mass — are considered to originate from the MS phase of these stars, 
have been observed in the optical (Lozinskaya 1982; Oey & Massey 1994; Marston 1997), 
in the IR (Marston 1996), and in the H I 21cm line (Cappa et al. 1996; Cappa & Benaglia 
1998; Benagha & Cappa 1999). 
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The MS shells are typical results of the interaction of massive stars with their sur- 
rounding ISM over a time period of a few Myr. There are also young objects which have 
influenced their interstellar environment so strongly that they can already be identified as 
star-gas interactions. One example that should certainly be mentioned here is the so-called 
Becklin-Neugebauer-Kleinman-Low (BN-KL) nebula in the Orion molecular cloud 1 (OMCl) 
behind the Orion nebula. It shows a spectacular, finger-like outflow of molecular hydrogen 
(Taylor et al. 1984; Allen & Burton 1993; McCaughrean & Mac Low 1997; Salas et al. 
1999). Due to extremely high visual extinction (Ay 20 — 50 mag) it can only be ob- 
served in the infrared or at longer wavelengths. Most models for the gas outflow proposed 
so far assume that it is driven by mass ejection from a young star, namely IRc2 and/or the 
Becklin-Ncugebauer (BN) object which are both very close to the projected center of the 
mass outflow. IRc2 has a mass-loss rate M„ ^ lO^*^^"'^) Mq yr^-*^ (Downes et al. 1981) and a 
luminosity Lphot ~ 2 ... 10 x 10'^ Lq (Genzel & Stutzki 1989). The mass-loss rate of BN is 
Mw ~ 4 X lO-^M0yr-i and the luminosity Lphot 1 . . . 2 x 10^ L© (Scoville et al. 1983). 

Another fairly well studied region where formation of stars and their interaction with 
the ambient medium takes place in an early phase is the Eagle Nebula (M16). Imaging with 
the WFPC2 onboard the Hubble Space Telescope (Hester et al. 1996) as well as observations 
in the infrared with ISOCAM (Pilbratt et al. 1998) and NIRSPEC (Levenson et al. 2000) 
reveal a number of interesting insights into the structure and possible formation history of 
this active region. The picture of M16 in the optical is dominated by the appearance of 
the so-called "elephant trunks", impressive pillars of dense molecular gas protruding into 
the H II region. Their surface is partly illuminated by the ionizing light of nearby (distance 
to the pillars ~ 2 pc) massive stars so that photo evaporative flows from the surface can be 
identified on the well resolved pictures (Hester et al. 1996) as striations extending normally 
from the surface of the cloud. The flow is driven by the pressure gradient that arises when 
the ultraviolet photons from the massive star heat the dense gas at the surface of the trunk 
to a few thousand Kelvin, resulting in a pressure which is higher than in the rest of the 
(lower density) H ll region. 

An interesting question arises in connection with M16: Is the elephant trunk structure 
a remnant of the original molecular cloud from which the young stars in the region were 
formed or has this structure actively been molded by the massive stars as a shell swept- up 
by stellar winds and folded by hydro dynamical instabilities? More detailed observations and 
better theoretical models are needed to answer this question. 
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2.4. Previous numerical models of stellar wind bubbles: an overview 

Since effects associated with inhomogeneous ambient media or time-dependent stellar 
parameters as well as the formation of gas dynamical instabilities are difficult to handle 
within an analytical framework, numerical simulations have become a powerful tool to study 
the evolution of H ll regions and SWBs. First 2D numerical calculations of SWBs have 
been performed by Rozyczka (1985) and Rozyczka & Tenorio-Tagle (1985a,b) in a series of 
papers. For the expansion of the bubble into a homogeneous medium the authors described 
the formation of clumps in the thin shell but due to the lack of resolution they were unable to 
completely resolve the dynamics of the shell and to compare it with theoretical predictions. 
The authors also studied how a change of the wind luminosity (an increase between two 
fixed levels in a certain time) infiuences the shell stability and how SWBs break out from 
a plane-parallel stratified disk. The shell fragmentation as result of an increase of the wind 
luminosity has also been studied by Stone, Xu, & Mundy (1995) and has been compared 
with the observation of "bullets" in the vicinity of young stars. 

GMLl and Garcia-Segura, Langer, & Mac Low (1996b) focused on the history of the 
stellar mass loss and the implications for the formation of WR bubbles. They calculated the 
evolution of the SWB in one dimension until the star enters the LBV phase (for the model 
with an initial stellar mass of 60 Mq) or until the star leaves the RSG stage (for the 35 M0 
model). The resulting ID profiles of circumstellar density, pressure and velocity were used 
to set up the ambient conditions for 2D calculations of the following stage when the WR 
wind interacts with the surrounding structure. The authors found that in both cases the 
slowly expanding shells originating from prior mass loss phases are heavily eroded when they 
are overtaken by the faster WR wind. This is also the phase in which the resulting nebula 
shell is considered to be most visible due to the high emission measure in the high density 
filaments resulting from the collision. The nebula NGC 6888 with its filamentary structure 
would fit in this picture as being driven by a WR star that has undergone MS and RSG 
evolution and the shell swept up by the WR wind is currently colliding with the RSG wind 
shell. 

The approach of Garcia-Segura & Franco (1996) is shghtly different because they were 
mainly concerned with the very early evolution of H ll regions. They presented 2D gas 
dynamical calculations for the evolution of H ll regions in constant and power-law density 
profiles and showed the effect of radiative cooling on the thickness of the swept-up shell 
and therefore also on the development of instabilities. They found that the ionization front 
can reinforce the growth of the thin-shell instability and in a power-law density fall-off with 
exponent equal to two it can even lead to violent shell disruption. 

Brighenti & D'Ercole (1997) investigated the formation of WR shells resulting from fast 
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WR winds evolving into slow anisotropic RSG winds. They found that the shell swept up 
by the WR wind becomes Rayleigh- Taylor unstable and fragments even before it hits the 
shell that had previously been piled up by the RSG wind. X-ray maps computed from the 
numerical models show that the anisotropy of the RSG wind leads to a "two-lobe" X-ray 
morphology which is qualitatively in agreement with the observation of NGC 6888 (Wrigge 
et al. 1994), but cannot reproduce the small filling factor of the X-ray emitting gas. 

Frank, Ryu, & Davidson (1998) tried a different approach: They modeled the interac- 
tion of an anisotropic fast wind with an isotropic slow wind to explain the morphology of 
LBV bubbles. They found that anisotropic fast winds can indeed produce strongly bipolar 
outflows without assuming that the fast wind collides with a slowly expanding disk or torus 
as previously postulated. 

Strickland & Stevens (1998) calculated simple MS bubbles blown by a wind with con- 
stant mass-loss rate and terminal velocity into a uniform ISM. Their main goal was to 
calculate synthetic X-ray spectra from these numerical models as they would be observed 
with the ROSAT satellite, in order to analyze these spectra according to standard procedures 
and to compare the inferred properties with the original numerical models. Surprisingly (or 
not), they found that the inferred properties of the bubble can considerably deviate from 
the "real" properties of the model bubble. This implies that detailed X-ray emission models 
are necessary — instead of simple one- or two-temperature spectral fits — in order to derive 
reliable properties of bubbles and superbubbles. 

In Table 1 we compare the numerical hydrodynamic models of SWBs / H ll regions 
quoted above with respect to the included physics, the covered range of parameters, and 
some technical properties. 

3. Numerical method 

3.1. The radiation- hydrodynamics scheme 

The numerical calculations presented in this paper have been performed using a 2D 
radiation-hydrodynamics code with a fairly sophisticated treatment of time- dependent ion- 
ization, radiation transfer, and heating and cooling processes. Most features of the code and 
results of test calculations have already been described in Yorke & Kaisig (1995) and Yorke 
& Welz (1996). Therefore, here we only give a short summary of the code's main properties. 

The hydrodynamical equations for the conservation of mass, momentum, and energy in 
an ideal, inviscid, single fiuid have been formulated in cylindrical coordinates {r,z), assuming 
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axial symmetry around the z axis. This set of equations is solved numerically on an Eulerian 
grid in the quadrant with r > and -2 > 0, i.e. we additionally assume mirror symmetry with 
respect to the equatorial plane. The differencing scheme used to discretize the equations is 
second order accurate in space. Due to operator splitting the accuracy in time is greater than 
first order. The advection scheme of van Leer (1977) is employed. The star that influences 
the gas is located in the center of the coordinate system at r = ;2 = 0. We use square 
grids with constant mesh size and insert multiply nested grids in the corner at r = z = 
to enhance the spatial resolution close to the star. All the nested grids are self-similar and 
the hnear spatial resolution is improved by a factor of two for each level of nesting. The 
basic code is exphcit, and therefore the Courant-Friedrichs-Lewy (CFL) condition determines 
the maximum time step except when changes in the stellar parameters take place on even 
smaller time scales. Von Neumann-Richtmyer artificial viscosity is used for the treatment 
of shocks. We use the equation of state for an ideal gas (in our case a mixture of molecular, 
neutral, and ionized hydrogen with electrons). The formation and dissociation of molecular 
hydrogen is not yet considered in the energy balance, because it is assumed that continuum 
radiation below the Lyman threshold dissociates the molecules before the gas is heated up 
to temperatures where thermal dissociation occurs. 

Our integration order starts with the finest grid. The sizes of the time step satisfies 
the CFL condition for the finest grid as well as for all coarser grids (normalized by the 
appropriate power of 2). After two time steps have been done on the finest grid and the 
solution has been advanced to ^ + 5^1 + 5^2, one time step Hx^bt^ is done on the next coarser 
grid to advance the solution there to the same time. Because of our prior careful choice of bt\ 
and (5^2, their sum obeys the CFL condition on the next coarser grid. After the solutions on 
the two grid levels have been advanced to the same point in time, all coarse grid values which 
have underlying fine grid values (not from the boundary) arc replaced by weighted averages 
from the fine grid. The outer boundary conditions of the fine grid are interpolated from the 
corresponding grid values of the coarse grid. After this cycle is repeated (taking into account 
the new CFL time-step constraints) the first integration on the next coarser (3rd level) grid 
can be performed; 3rd level values are partly replaced by 2nd level values and 2nd level outer 
boundaries are imposed from 3rd level values, and so on. We use this recursive scheme for 7 
grid levels, which means that X^^^q ~ ^ ~ 1) — 127 single integrations have to be done 
before one time step on the coarsest grid is completed. 

The ionization structure of hydrogen is calculated for each hydrodynamical time step 
by including the effects of photoionization, coUisional ionization, and spontaneous recombi- 
nation. The absorption of the stellar point-source Lyman continuum by neutral hydrogen 
is calculated along radial lines-of- sight. We use iVr -f- iV^ — 1 rays to ensure that every grid 
cell is traversed by at least one ray, where A^r and A^ are the numbers of grid cells in r and 
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z direction, respectively, without ghost cells. The photoionizing radiation field is assumed 
to have a black body spectrum at the effective temperature of the star. We use the grey 
approximation for the transfer of the stellar photons, i.e. we calculate the absorption cross 
section and excess energy only for the mean energy of photons above the Lyman threshold 
and, therefore, do not account for selective absorption processes. To account for the dif- 
fuse Lyman continuum that originates from electron recombinations directly into the ground 
state of hydrogen we simply use the on-the-spot approximation. 

Radiative heating and cooling are considered as source and sink terms in the en- 
ergy equation and arc calculated in an additional substcp. Photoionization of hydrogen 
is responsible for gas heating. The treatment of cooling depends on the gas tempera- 
ture. Below 15 000 K the contributions by the following processes are exphcitly calculated: 
Bremsstrahlung, coUisional ionization of hydrogen, thermal energy lost in hydrogen recom- 
binations, coUisional excitation of neutral hydrogen followed by Lya emission, coUisional 
excitation of the low-lying terms of N"*" and O"*"^ and the term of O"^. Simple as- 
sumptions have been made about the degree of ionization of nitrogen and oxygen: Due to the 
similar ionization potentials and comparable recombination coefficients of H*^, N*^, and 0°, 
we set the ratios of ionized-to-neutral species for nitrogen and oxygen equal to the respective 
number calculated for hydrogen, i.e. Pn+/pn = ^ and (po+ + Po++)/Po = where X is 
the hydrogen ionization fraction. It is furthermore assumed that the oxygen ions are equally 
distributed on the ionization stages O"^ and O"^"^, i.e. po+ 1 Po = Po++ 1 Po = X/2. 

The chemical composition is solar everywhere in our computational domain. Helium 
is not yet considered because it is much less abundant than hydrogen and has no low-lying 
energy levels. For temperatures above 10^ K we use the interstellar cooling function from 
Sarazin & White (1987) with the correction given by Soker (1990). This cooling function is 
based on the assumption of coUisional ionization equilibrium. In the intermediate tempera- 
ture range between 15 000 K and 10^ K we calculate both values and use a weighted average. 
We have not included dust in the calculations presented here. 

3.2. Initial conditions 

We start all our calculations with the turn-on of the ZAMS stellar radiation field and 
stellar wind in a homogeneous and quiescent ambient medium. This is definitely a gross 
oversimplification of what can be expected for the structure of the circumstellar gas shortly 
after the star has been formed at the end of the prc-MS phase. We have two reasons for 
our approach: The first is that, in order to understand what happens in the circumstellar 
medium, it is important to start with a well defined and simple initial configuration. As 



- 12 - 



we will show below, even with this initially homogeneous ambient density and temperature 
distribution, the interaction of the SWB with the ionizing radiation field produces a variety 
of interesting morphological structures whose formation processes have to be understood 
in detail before studies with a more realistic initial setup can be performed. The second 
reason is that wc want to keep our simulations comparable with those of other authors who 
use different algorithms or incorporate different physical effects in their codes. We choose 
no — 20 cm~^ and Tq = 200 K for our quiescent ambient medium. The density is perhaps high 
for a star that has already left its parental molecular cloud, but we want to keep our results 
comparable to those of GMLl who also used uq — 20cm~^ in order to reduce computational 
expenses. The temperature Tq = 200 K in the ambient medium was chosen to yield a thermal 
pressure "typical" for the ISM at the solar galacto-centric distance. 

We use a spherical ID solution after 1000 yr as initial model for the 2D simulations 
in order to prevent boundary effects from influencing the formation of the very small wind 
bubble while the development of the self-similar structure is not completed. We switch from 
ID to 2D before the swept-up gas collapses into a thin shell that is subsequently subject to 
instabilities. 



3.3. Boundary conditions 

In our models the stellar parameters for mass-loss rate (M^), terminal velocity of the 
wind (fw), effective temperature (Tes), and photon luminosity in the Lyman continuum 
(LLyc) are time-dependent boundary conditions which drive and govern the evolution of 
the circumstellar gas. For the 60 star we adopt the stellar parameters given by GMLl. 
They were obtained from stellar evolutionary models and from observations and are shown 
in Figure 2. The 60 M© star is supposed to undergo the following evolution: MS O star —>■ 
H-rich WN star P Cygni-type LBV H-poor WN star H-free WN star WC star 
^ SN. 

According to the mass-loss rate and terminal velocity at time t, appropriate values for 
density and velocity of the gas are set within the "wind generator region" on the finest grid 
(a small sphere around the center of the coordinate system where the star is located). The 
radius of the "wind generator region" is 3.5 x 10^^ cm for the calculations shown in this 
paper. 

Reflecting boundary conditions are used at the r- and 2;-axis. The outer boundaries of 
the nested grids were taken from the corresponding values on the next coarser grid. The 
outer boundaries of the outermost (coarsest) grid are semi-permeable, i.e. mass, energy. 
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and momentum can flow out of the computational domain, but not in. Of course, outflow 
would conflict with our intention to carefully take stock of the energetic processes in our 
calculations. Thus, we prevent outflow by choosing the coarsest grid large enough that 
neither the ionization front nor moving material can reach the outer boundary before the 
star dies and the calculation stops. 

3.4. Geometry and resolution 

The size of the computational domain is rmax = -^max = 64 pc. We use 6 nested grids 
within the coarsest grid, resulting in 7 grid levels. The simulation has been performed with 
125 X 125 cells on each grid level (excluding ghost cells). The linear resolution ranges from 
8 X 10^^ pc close to the star to approximately 0.5 pc in the outermost parts of the coarsest 
grid which were affected only during the late stages of the bubble evolution. For a resolution 
study we repeat the calculations with 61 x 61 cells and for the first Myr of evolution with 
253 X 253 cells. 



4. Results and discussion 

4.1. A resolution study 

Before we present our results we briefly discuss their validity. It is currently impossible 
to spatially resolve all details of the entire dynamical evolution when modeling an SWB 

in two or more dimensions over such a long time with all the physics we have included. 
Especially the modes of the thin-shell instability (Vishniac 1983; Ryu & Vishniac 1988; 
Vishniac & Ryu 1989) arc extremely difficult to resolve, because a sufficient number of grid 
cells across the thin shell is required to follow the tangential flow of material. Mac Low 
& Norman (1993) performed purely hydrodynamical simulations of the thin-shell instability 
and confirmed the linear stability analysis of Vishniac (1983). Since they completely focused 
on the hydrodynamical evolution of the thin shell, they made no efforts to treat much more 
complex systems with more sophisticated physics. 

To check the impact of this resolution effect on our results for the energetic evolution of 
the SWB, we performed a resolution study. We calculated the same 60 M0 model with three 
different resolutions, medium (125 cells in each dimension), low (61 cells in each dimension), 
and high (253 cells in each dimension). We compare only the first Myr of the evolution 
(because the high-resolution model is too CPU intensive). The results are displayed in 
Figure 3. The change of kinetic energy due to different resolutions is systematic (at times 
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later than 0.2 Myr) in the sense that a higher resolution increases the kinetic energy that 
remains in the system. For most of the time the difference in kinetic energy between medium 
and high resolution (< 0.05 dex) is smaller than between low and medium resolution, which 
is expected as the value converges toward the limit. The variation of the thermal energy with 
resolution is only weak, the maximum difference between low and high resolution is < 0.05 
dex at any time during the first Myr. The scatter in ionization energy for the three different 
resolutions is also < 0.05 dex for 0.5 Myr < t < 1.0 Myr and < 0.1 dex for t < 0.5 Myr. In 
the latter case, the deviation of ionization energy between the medium and high resolution 
runs varies, but is considerably smaller than between low and medium resolution, indicating 
that the value of the ionization energy is already close to the actual value. 

These results show that the errors in our energetic analysis due to resolution effects are 
within an acceptable range, although some details of the evolution are not yet completely 
resolved in our calculations. 

4.2. The evolution of the 60 Mq model 

Figures 4 to 16 show the evolution of the circumstellar gas around our 60 Mq model 
star. The plots show grey-scale coded physical quantities in a quadrant of the r-z-plane 
through the star. The size of the displayed area varies with time as the expansion proceeds 
outward. The star is located in the center of the coordinate system. We display the mass 
density together with the velocity field, the gas temperature, and the degree of hydrogen 
ionization, respectively, for certain stages of evolution. The degree of ionization is not shown 
for stages after 1 Myr, when the ionized region spatially almost coincides with the hot gas. 
To prevent confusion, velocity arrows are skipped in the free-flowing wind and (except for 
Figure 4) in the hot bubble. 

Figure 4 shows the initial model for the 2D calculation that has been set up from the ID 
solution after 1000 yr. The fundamental structure of the combined SWB / H II region is well 
described by the analytical solutions. The ionization front at a distance of about ~ 8pc from 
the star separates the ionized gas at about 8000 K from the undisturbed medium. The front 
has not yet reached the initial Stromgren radius at about 12.7 pc and is still R-type. Thus, 
there is httle dynamical response of the heated gas at this point. The SWB evolves into 
the pre-ionized and pre-heated but otherwise undisturbed medium. The free-flowing wind is 
heated by the reverse shock at r ^ 0.16 pc to very high temperatures of about 10^ K. Though 
the density in the hot gas is fairly low (fsi 4 x 10~^^gcm~^), the thermal pressure is still 
much higher than in the H II region and the hot bubble expands supersonically (with respect 
to the isothermal sound speed in the H ll region of ~ 12kms~^) at roughly 200kms~^. The 
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photoionized H ll region gas becomes shocked by the expanding hot bubble and compressed 
into a shell with a density of more than 10~^^ gcm~^ and a temperature of nearly 10^ K. At 
this early time, the SWB and the photoionized H ll region still evolve independently. 

Figure 5 depicts the evolution after 2 x 10^ yr. The photoionized region has reached 
its initial radiative equilibrium size and gas at the outer edge of the H ll region starts to 
accelerate outward due to the pressure gradient between the warm ionized gas and the cold 
neutral ambient medium. The hot bubble has grown to about 2 pc in radius. The shell has 
been decelerated to approximately 60kms~^. Due to the high density and high temperature, 
strong cooling has led to compression of the swcpt-up gas into a thin shell. The thin-shcll 
instability initially produces some dense knots in the shell that imply variations in the optical 
depth along radial lines-of-sight indicated by the rippling of the photoionization front. The 
initial evolution of the SWB (as long as it is independent of the H ll region evolution) is 
very similar to what is described by Strickland & Stevens (1998), except for the deviations 
due to the different ambient density and stellar parameters. 

On the basis of Figure 6 {t — A x 10^ yr) and Figure 7 {t — 0.1 Myr) the structuring 
impact of the interaction processes between the SWB and the H ll region can be well iden- 
tified: The photoionized H ll region has started to expand into the ambient medium at 6 to 
7kms~^. The isothermal sound speed of the cold neutral gas is about lkms~^. Thus, an 
outer radiative forward shock occurs that sweeps up the ambient gas into a shell. The hot 
gas bubble has grown to about 5pc in radius and the shell expansion velocity is still about 
25kms~^ (t = 0.1 Myr). The ionization plots nicely illustrate the influence of the density 
knots in the wind bubble shell on the evolution of the H ll region: The gas clumps cast 
shadows into the H ll region. The optical depth for hydrogen ionizing photons along lines- 
of-sight through the clumps is higher than for lines-of-sight which pass the shell between 
clumps. Due to the weakening of the radiation field, the gas behind the clumps starts to 
recombine, forming neutral spokes within the H ll region. The spokes are less prominent in 
temperature than in degree of ionization, because the cooling efficiency breaks down with 
the disappearance of the free electrons and the shadowed regions remain at several thousand 
K for some time. 

One could imagine that the shadowed regions are artificial due to the use of the on-the- 
spot approximation for the diffuse radiation field that originates from recombinations directly 
into the ground state of hydrogen. But a test calculation that used fiux limited diffusion 
instead of the on-the-spot approximation yielded qualitatively the same results with respect 
to the shadowed regions. Only slight morphological differences were obtained due to the 
transfer of the diffuse radiation field. 

Considerations about the infiuence of shadow patterns on ionized regions are not com- 
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pletely new: Capriotti (1973) already discussed the impact of shadows by shell fragments on 
the formation and appearance of planetary nebulae (PNe). Williams (1999) examined the 
corrugation of R-type ionization fronts by density inhomogeneities with subsequent forma- 
tion of dense clumps. Soker (1998) studied the formation of compressed tails in the shadow 
zones behind dense clumps which are present in the vicinity of the central star of a PN before 
the ionization starts. And Canto et al. (1998) presented an analytical model for the partial 
ionization of the shadow regions behind neutral clumps by the diffuse radiation produced in 
the nebula along with some gas dynamic simulations. The model is not restricted to PNe 
and can be generalized to photoionized regions. 

Richling & Yorke (2000) considered the effects of UV dust scattering and calculated 
the evolution of disks externally illuminated by an O star. Here, also, shadows are cast 
which appear cometary in shape with tails pointing away from the ionizing source. Had 
we considered UV scattering by dust in this investigation it is very likely that our clumps' 
shadows would also have had a cometary shape. 

For the combined H ll region / SWB described in our work we show that the density 
inhomogeneities needed to cast the shadows can be produced by the action of the stellar 
wind even after the H ll region has already been formed in a completely uniform medium. 

In our case the growth of the ionized fingers of the H II region is therefore not only due to 
the ionization front instability described by Garcfa-Segura & Franco (1996). It is strongly 
triggered by the redistribution of mass (and thus opacity) by the action of the stellar wind 
shell. 

The real nature of the spokes, clumps, and instabilities is 3D rather than 2D, a fact which 
certainly should infiuence their size distribution, shapes, and longevity. 3D simulations, 
though extremely expensive, would be useful to study these structures in more detail. These 
2D simulations should thus be regarded as qualitative indicator that such instabilities exist, 
a feature missing from ID calculations. 

Another consequence for the structure of the H II region and the shell of the SWB can 
be seen on Figure 8 at 0.22 Myr. The disturbed structure of the H II region, especially the 
formation of the neutral spokes and the ionized fingers, has induced several fiow patterns with 
velocities up to ^ 20 km s~^. Due to the pressure gradient between the neutral spokes and the 
ionized gas, the spokes are compressed by the ionized gas. At the same time the expanding 
SWB compresses the H ll region (and its spokes) as a whole. The clumps in the stellar 
wind shell originally responsible for the shadows grow in mass and partially merge during 
the further coiusc of evolution. Photoevaporation of the stellar wind shell and its embedded 
clumps docs not strongly infiuence the dynamics, because the gas in these structures is still 
warm and its photoionization does not produce large pressure gradients. However, where 
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stellar photons directly impinges upon the cold outer shell through the optically thinner 
ionized fingers photocvaporation occurs and leads to significant dynamical effects, as evident 
in Figures 8 through 16. 

In the top panel of Figure 8 and also in some of the other density plots are swirls of 
material visible in the hot bubble. These vortices are sometimes generated at the r- and 
2;-axes, possibly as an artifact of the refiecting boundary conditions. Since the total mass 
involved in these flow patterns is very small, they do not influence the morphological or the 
energetic evolution signiflcantly. 

Proceeding to 0.5 Myr (Fig. 9, note the different scale) we see that the entire H ll 
region has been swept up by the SWB, which has grown to approximately 12 pc in radius. 
The photoionized H II region (it can best be distinguished from the hot bubble using the 
temperature plot) consists of 3-7 pc long pillars, the remnants of the ionized fingers. The 
stellar wind shell loses its structure, becomes more or less dissolved, and merges into the 
"remnant" H ll region, because it sweeps up a highly nonuniform medium with an expansion 
speed that becomes comparable to the sound speed of the H ll gas. 

Figures 10 and 11 compare the results of our high and medium resolution runs with 
respect to the morphological structure after 1 Myr. One can easily see that there are signif- 
icant deviations in the structure of the photoionized region between the two runs. This is 
not very surprising, because, as mentioned above, the formation of density fiuctuations due 
to the thin-shell instability is very sensitive to the spatial resolution, especially as we have 
started with a completely undisturbed ambient medium and have not artificially excited any 
particular unstable mode. And we have also seen that exactly these clumps in the stellar 
wind shell strongly influence the evolution of the H ll region due to their opacity for Lyman 
continuum radiation. 

Nevertheless, the general structure is very similar for both resolutions: The radius of the 
hot bubble is roughly 17 pc (strong deviations from spherical shape appear in both cases). 
The photoionized gas forms a highly irregular shell around the hot bubble with a mean 
thickness of ~ 5 pc, bound on the outside by a swept-up shell of shocked ambient material. 
It appears as if the phenomenon of "break-out" of the hot gas into the photoionized gas has 
occurred between the spokes in both cases. However, we see no evidence for further break- 
out into the undisturbed constant density medium. Had we chosen an ambient environment 
of limited size — i.e. included cloud edges of decreasing density within the computational 
grid — we could expect a more dramatic demonstration of break-out. 

We now jump to t — 3.30 Myr and show in Figure 12 the evolutionary state of the cir- 
cumstellar material shortly before the star enters the LBV phase with its strongly enhanced 
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mass-loss rate. (Please note that we once again changed the plot scale to an appropriate 
size.) The radius of the hot bubble has grown to approximately 40 pc. Due to the growth 
of the bubble (and the slight decrease of H-ionizing photon luminosity after 2 Myr) the pho- 
toionized H II region has more or less collapsed into a skin with ^ 2 pc thickness at the inside 
of the outer shell. The outer shell itself is still corrugated, and so is the photoionized gas. 
Some small clumps have become detached from the shell. 

Figure 13 depicts the structure during the LBV eruption. Compared with the earlier 
plots one can detect a much higher density in the free-flowing wind zone. That is due to the 
higher mass-loss rate of the star and the smaller terminal velocity. The rest of the structure 
is largely unchanged, because the time-step difference is small. 

The LBV phase lasts only for about 10 000 yr and is followed by the final WR phases of 
the star, which are once again associated with wind terminal velocities of a few 1000 km s~^. 
GMLl describes what happens next: The WR wind becomes shocked when it hits the much 
slower LBV wind and the hot, rarefied gas accelerates the dense LBV wind (Fig. 14). This is 
now a highly Rayleigh- Taylor unstable situation and the shocked WR wind breaks through 
the LBV material, fragments it, and blows it into the MS bubble (Fig. 15) where it mixes 
with the hot, rarefied gas, enhancing the density and lowering the temperature there for 
some time. 

The last set of pictures (Fig. 16) shows the structure of the circumstellar gas at the end 
of our simulation. This is the pre-SN configuration immediately before the star explodes 
as SN II. The radius of the hot bubble has reached nearly 50 pc and the photoionized H ll 
region is limited to the illuminated inner skin of the outer shell of swept-up ambient gas. 

A "re- ionization" phase (cf. Beltrametti, Tenorio-Tagle, & Yorke 1982; Tenorio-Tagle et 
al. 1982) cannot be found in our model, because the dynamical evolution of the circumstellar 
gas at late stages is dominated by the infiuence of the stellar wind rather than the stellar 
radiation field. 



4.3. The energy balance in the circumstellar gas 

4-3.1. Numerical results 

Figure 17 shows the total energy contributions in the circumstellar gas of the 60 Mq 
model as a function of time. We follow the kinetic energy of bulk motion throughout the 
computational domain, the ionization energy (13.6 eV per ionized hydrogen atom), the ther- 
mal energy of cold (T < lO^K), warm (lO^K < T < 10^ K), and hot (T > 10^ K) gas. As 
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we have seen before, the ionization and heating of the H ll region is the fastest process at 
the beginning of the evolution. Thus, the ionization energy reaches 1.1 x 10^° ergs and the 
thermal energy of warm gas around 1.5 x 10^^ ergs very soon after the stellar turn-on. 

The evolution of the thermal energy of warm gas follows that of the ionization energy 
over the lifetime of the star with a shift of 0.7-0.9 dex. This is due to the fact that photoion- 
ization by the stellar radiation field is responsible for the bulk production of ionized gas at 
typically 8000 K. The contribution by shock heating is small, because the reverse shock does 
not heat much mass and the thermal energy of this gas is accounted under "hot gas" rather 
than "warm gas". Therefore, we will discuss only the ionization energy bearing in mind that 
the thermal energy of warm gas behaves quite similarly. 

An interesting feature is the drop of ionization energy between its first maximum and 
t ~ 0.2 Myr. This is due to the formation of the dense structures (shell, spokes) in the 
H II region which reduce the recombination time and increase the cooling of the gas (see 
also Fig. 18). Between t 0.2 Myr and t ^ 1.9 Myr the denser structures dissolve and 
the pressure decay within the hot bubble together with the H ll region expansion result in 
a decrease of density in the H ll region. Thus, the recombination time increases and the 
ionization energy stored in the system increases accordingly. The ionization energy reaches 
2.7 X 10^° ergs and the thermal energy of warm gas 3.7 x 10^^ ergs. 

As the star enters the first WN stage at i 1.7 Myr, the photon luminosity in the Lyman 
continuum begins to decrease slightly while the mechanical wind luminosity starts to increase 
(see Fig. 2). The former directly reduces the rate of photons available for photoionization, 
whereas the latter enhances the pressure in the hot bubble, thus compressing the H ll region, 
which implies higher recombination rates. Thus, ionization energy and thermal energy of 
warm gas decrease by a factor of 2 and reach local minima at t ^ 2.6 Myr. They rise again 
when the photon luminosity in the Lyman continuum increases and the mechanical wind 
luminosity decreases until the star reaches the LBV stage. 

Though spectacular in appearance, the impact of the LBV phase on the energy balance 
is only limited due to a relatively brief period of time. We note here that the table of stellar 
parameters that we use is sufficiently time resolved even during the LBV phase, contrary to 
what could be suspected from Figure 2. For a detailed plot of the stellar parameters during 
the LBV phase see GMLl. The ionization energy in the circumstellar gas drops during 
the LBV stage, because the effective temperature of the star decreases significantly. But 
immediately after the LBV eruption during the next WN stage the stellar UV flux drives 
the circumstellar ionization energy to a global maximum at 3.4 x 10™ ergs. The ionization 
energy follows the stellar UV flux. The flnal values of ionization energy and thermal energy 
of warm gas (i.e. at the end of our calculation) are 1.0 x 10^° ergs and 2.1 x 10"^^ ergs. 
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respectively. 

The kinetic energy of bulk motion in the circumstellar gas rises from zero at the begin- 
ning (because we started with a quiescent medium) as more and more circumstellar gas is 
accelerated by the expansion of the H ll region and the SWB. After some 0.3 Myr it reaches 
the same level as the thermal energy of warm gas and rises in the same manner. When the 
mechanical luminosity of the stellar wind increases in the first WR phase it lifts the kinetic 
energy above 10^'^ ergs, close to the ionization energy level. The ejection of the LBV nebula 
raises the kinetic energy for a short time almost by a factor of 2, but it drops back when the 
LBV material hits the outer shell and the kinetic energy is dissipated into thermal energy. 
Its final value is lA x 10^'' ergs. 

The thermal energy of hot gas accounts for the internal energy of the rarefied gas in 
the SWB that is heated to high temperatures by the reverse shock during the MS and WR 
phases of the star. Since there is no hot gas at the very beginning, the thermal energy of 
this gas phase rises from 0. The behavior of the curve is very similar to that of the kinetic 
energy but it increases faster during the first WN stage of the star. Between 2.3 and 3.1 Myr 
it is the dominant form of energy in our system and after some slight changes in and around 
the LBV phase it finally decays with the mechanical luminosity of the star to 9.7 x 10^^ ergs. 

The thermal energy of cold gas starts already with the internal energy of the quiescent 
ambient medium. The value increases smoothly from 3.4 x 10^^ ergs at the beginning to 
7.9 X 10^^ ergs at the end of the calculation, because more and more ambient gas becomes 
swept up by the outer shell. The weak shock heats it to several 100 K whereas coohng of 
this neutral gas is very low. 

At the beginning the radiative energy input has the strongest infiuence on the energy 
distribution. The dynamical response of the circumstellar gas (the expansion of the H ll 
region and the acceleration of the material by the SWB) takes much longer. Thus, the 
ionization energy dominates the circumstellar energy distribution over the first ~ 2 Myr. 
After the kinetic energy and thermal energy of hot gas have caught up with the thermal 
energy of warm gas, these three forms of energy are comparable between 0.3 and 1.8 Myr. 
Later, due to the accumulation of kinetic energy in the shell and thermal energy of hot gas 
in the bubble, these two become comparable with the ionization energy. Between 2.2 Myr 
and the end of the calculation these three forms of energy differ by not more than a factor 
of 3, though some changes occur. The LBV phase of the star induces very rapid changes in 
the relative energy ratios of at most a factor of 3 but they are transitory. 

It is interesting to compare the energy stored in the system for the two cases: 1) with 
wind and 2) the pure H ll region evolution without wind (Fig. 18). Obviously, the total kinetic 
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energy in the calculation with wind is always higher than in the corresponding calculation 
without, because of the added kinetic energy of the SWB shell. At the end of the calculation, 
the kinetic energy in the model with stellar wind is 4 times higher than in the model without. 

By contrast, the ionization energy reaches the same level oi ^ 1.1 x 10^° ergs very soon 
after the start of the calculation in both cases. As we have already seen in Figure 17, the 
ionization energy in the model with wind drops by a factor of 2 during the first 0.2 Myr 
whereas in the model without wind it continues to grow. This supports our explanation that 
inhomogeneities in the H II region (shell, spokes) result in shorter recombination times and 
thus lead to fewer photoionized hydrogen atoms in the system. Although there is additional 
energy input by the stellar wind, the ionization energy is lower than in the calculation without 
wind by 0.1 to 0.4 dex for most of the evolution. This, however, does not imply that the Ha 
luminosity is different. The Ha luminosity directly mirrors the stellar ionization input. 

Although the thermal energy of warm gas (55 % below the respective value of the wind- 
less model at the end of the calculation) is similarly affected as the ionization energy, the 
total thermal energy is always higher in the model with wind, because we have additional 
energy deposited into the hot and into the cold gas phases. 

Whereas the mechanical power of the stellar wind and the radiative power of the star 
are the energy sources in our system, the only energy sink term is radiative cooling of the 
gas. All other energetic processes lead to a redistribution of forms of energy within the 
system but not to a change of the total energy. Figure 19 shows the total energy loss of the 
system due to cooling of the gas and distinguishes 3 different contributions: "H ionization 
loss" refers to the 13.6 eV binding energy per hydrogen recombination into all states above 
the ground state, while "H thermal loss" accounts for the thermal energy of the electrons 
lost by the same process. "Other processes" include all other cooling processes considered 
such as collisionally excited line emission and Bremsstrahlung. One can see from Figure 
19 that the shape of all the curves follows the Lyman continuum photon luminosity of the 
star (sec Fig. 2). The total energy-loss rate during the MS phase is close to 10^^ ergs s~-^, 
approximately the same as the stellar Lyman continuum luminosity. Although the ionization 
energy in the system rises to 10^° ergs during the first few 10 000 yr, it is only a small fraction 
of the total input energy. 

"H ionization loss" is the dominant energy sink term during the evolution. It accounts 
for roughly 80 % of the total emitted power during the MS and first WN stage of the star. 
After the LBV phase, the relative importance of "other processes" has risen because the 
LBV ejecta enhanced the density in the hot bubble where collisionally excited lines and 
Bremsstrahlung dominate the energy emission. 
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Figure 20 shows the transfer efficiency into kinetic, ionization, and thermal energy and 
their sum for the model without wind. We define the transfer efficiency as the cumulative 
fraction of the input energy that has been converted into a particular form up to the time 
t — T. In this case we have only radiative input energy and all values are related to that. 

The transfer efficiency into ionization energy has the highest value during the whole 
calculation. The transfer efficiencies into ionization and thermal energy strongly drop at 
the beginning of the calculation. This is a consequence of the fact that the photoionization 
equilibrium has been quickly established and the ionization energy stays constant within 
an order of magnitude, whereas the transfer efficiency is related to the accumulated input 
energy that grows continuously in time. The transfer efficiency into ionization energy reaches 
5 X 10~^ immediately before the LBV phase and drops due to the decreasing photon 
luminosity in the final WR phases by almost a factor of 5. The transfer efficiency into 
thermal energy reaches the level of 10~^ before the LBV phase, but decreases only by less 
than a factor of 2 during the remaining evolution. The discrepancy between the declines of 
both values can again be explained by the fact that cooling in formerly photoionized gas is 
greatly reduced due to missing electrons after recombination occurs. Thus, the recombined 
neutral gas has lost all its ionization energy but can retain a considerable portion of its 
thermal energy over a relatively long period. The transfer efficiency into kinetic energy has 
reached a nearly constant level of 3 x 10~^ after 1 Myr and remains at that level until the 
end of the calculation, i.e. after the first Myr the kinetic energy in the system grows almost 
proportionally to the total input energy. 

Now we check the impact of the stellar wind on the energy transfer in the system. 
Figure 21 shows the same transfer efficiencies as Figure 20 but for the model with wind 
and radiation. Whereas for the model without wind the transfer efficiency into ionization 
energy dominated the total energy transfer efficiency, for the model with wind this is only 
true during the first ~ 2.2 Myr. In the first WN phase the transfer efficiency into thermal 
energy dominates for about 1 Myr, because the photon luminosity in the Lyman continuum 
decreases and the mechanical wind luminosity increases, thus enhancing the production of 
hot gas. Curiously, in the LBV phase all three energy transfer efficiencies have the same 
value of about 2 x 10~^. As in the case without wind, the transfer efficiency into ionization 
energy drops more strongly after the LBV phase than that into thermal energy. The drop of 
the latter one in this model is mostly due to energy loss of the hot gas that dominates the 
thermal energy at late stages of the evolution. 

The total energy transfer efficiency reaches approximately 4 x 10"'^ at the end of the 
calculation, about twice as much as in the model without wind. This is interesting, because 
we know from Figure 2 that the mechanical wind luminosity integrated over the whole 
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lifetime of the star is almost negligible compared to the total radiative energy input in the 
Lyman continuum. In other words, although the stellar wind itself does not actually inject a 
considerable amount of energy into the circumstellar gas (compared to the stellar radiation 
field), its presence almost doubles the total energy that is finally contained in the gas. This 
difference is due to the fact that the thermal energy deposited increased by a factor of 2.4 
and that the kinetic energy rose by a factor of 4 compared to the model without wind. The 
ionization energy actually decreased shghtly by some 20 — 30 %. 

To facilitate the discussion in the next section we summarize the values of the individual 
energy components at the end of the simulation in Table 2. Besides the kinetic energy of 
bulk motion (E'k) and the ionization energy of hydrogen {E[) we list the thermal energy of 
the cold, warm, and hot gas components (-E't,coid) -E't,warm, and -Et,hot) for the two models, one 
with and one without stellar wind. The values of the energy transfer efficiency into kinetic 
energy (sk), ionization energy (si), and thermal energy (et) at the end of the simulation are 
given in Table 3. The total energy of the Lyman continuum radiation emitted by the star 
is -ELyC = 1-07 X 10^^ ergs and the mechanical energy injected into the system by the stellar 
wind, if considered, amounts to — 3.31 x 10^^ ergs. 

4-3.2. Comparison with analytical results 

We now calculate the values for the kinetic-, ionization- and thermal energy in the 
system according to the analytical solutions given in section 2. The analytical approach 
cannot handle time-dependent stellar parameters, therefore we simply choose mean values 
for effective temperature and luminosity in the Lyman continuum over the lifetime of our 
model star. With (T^s) = 5.03 x 10^ K, (LLyc) = 8.33 x 103^ergss-\ and using = 
3.37 X 10~^^cm^s~^ as hydrogen recombination coefficient and Cs,n = 1.15 x lO^cms"^ for 
the isothermal sound speed in the H ll region (corresponding to T„ = 8000 K), we obtain 
for no — 20cm~^ after r = 4.065 Myr from equations (3), (4), and (5) for the model without 
wind: 

= 4.3 X 10^^ ergs , 
Ei = 5.0 X 10^° ergs , 
Et = 7.7 X 10^^ ergs . 

The energy transfer efficiency in the analytical approach is then defined as 

E 

£ — 



(11) 
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where E can be any of E]^, Ei, and Et, depending on the efficiency that is calculated. Thus, 
we obtain 

£k = 4.0 X 10~^ , 
£; = 4.7 X 10"^ , 
£t = 7.2 X 10"^ . 

Bearing in mind all the assumptions and approximations which have been made to 
obtain the analytical expression for the kinetic energy of bulk motion (see Lasker 1967), the 
deviation of less than 30 % in £k is not too bad. This difference is not due to the constant 
effective temperature and luminosity that we have chosen to obtain the analytical result. 
Recalculating for the first Myr, where the stellar parameters are about constant, makes 
the discrepancy between the analytical and numerical value of even larger. The difference 
is more likely due to temperature deviations from 8000 K in the H ll region of the model 
calculation. 

Comparing the analytical and numerical results for the transfer efficiency into ionization 
energy shows fairly good correspondence immediately before the LBV phase. At the end of 
the calculation the numerically determined transfer efficiency is about 0.6 dex below the an- 
alytical value because the Lyman continuum luminosity of the star drops significantly during 
the final WR phases of the star and the ionization energy follows immediately. Of course, 
the analytical solution cannot reproduce this feature. The correlation between the drop 
of Lyman continuum luminosity and the thermal energy is weaker, because cooling ceases 
when the plasma becomes neutral, as already mentioned above. Thus, the correspondence 
between analytical and numerical transfer efficiency into thermal energy is much better and 
the deviation at the end of the calculation is less than 15%. 

It is more difficult to compare analytical with numerical results in the case of the com- 
bined H II region / SWB model because we have only the analytical energy transfer efficiency 
solutions for the H II region and SWB separately. We simply add up the energy contributions 
from the H II region and the SWB bearing in mind that this is only a rough approxima- 
tion, which actually neglects the interactions between both structures. Because we have not 
considered cooling in the hot bubble, the analytical energy transfer rates into kinetic and 
thermal energy are upper limits. 

We insert the mean value of the mechanical wind luminosity, (L^) = 2.58 x 10^^ ergs s^^, 
into equations (9) and (10) and together with the results for the pure H ll region we obtain 

El, = 9.5 X 10™ ergs , 
Ei = 5.0 X 10™ ergs , 
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Et = 1.6 X 10^' ergs . 

To find the energy transfer efficiency, we relate tfiese values to the sum of Lyman continuum 
radiation energy and the mechanical wind energy (which, of course, is almost negligible), 

_ E 

' ~ r iiUyc) + (K)) ' ^^^^ 
where E can again be any of £^k, -^i, and E^. We get 

£k = 8.6 X 10-^ , 
£i = 4.6 X 10"^ , 
£t = 1.4 X 10"^ . 

We have already discussed the fact that the ionization energy in the numerical model with 
wind is lower than in the model without wind. This slightly enhances the discrepancy 
between the analytical and the numerical results for the ionization energy in the case with 
stellar wind. On the other hand, there is no indication that the neglect of the ionization 
energy in the hot bubble is a bad approximation. 

The increase of the thermal energy deposit from the windless H ll region model to the 
combined SWB/H ll region model (factor 2.4) is almost an order of magnitude below the 
analytical upper limit. Since the conversion of thermal energy by PdV work is the major 
source of kinetic energy, the increase of the kinetic energy deposit between the models without 
and with wind (factor of 4 at the end of the calculation) is also 82 % below the analytical 
upper limit. These findings arc supported by the observation that the radius of the bubble 
at the end of the calculation 50 pc) is considerably smaller than the analytical result of 
~ 68 pc according to equation (6). One reason for the lack of thermal energy in the numerical 
model is the variation of the mechanical luminosity of the star. The LBV wind enhances the 
density in the bubble which leads to stronger cooling, and the mechanical luminosity of the 
stellar wind generally decreases during the last 0.6 Myr. But also during the other stages 
of evolution, coohng in the hot bubble is not completely negligible for the energetics of the 
system. The thermal energy of the hot bubble calculated analytically t — 1 Myr is already 
a factor of 4.5 larger than the corresponding value in the numerical model. Although we have 
not implemented heat conduction, effects like turbulent mixing between hot and cold gas as 
well as numerical diffusion enhance the cooling rate. The resolution dependent numerical 
diffusion is not a key factor as wc have shown in our resolution study. However, the cooling 
due to gas mixing is enhanced because the SWB expands at the beginning into the highly 
nonuniform H II region, a fact that cannot be properly handled by the analytical description. 
In general, the analytical results for the upper limits of the transfer efficiency into thermal 
and thus kinetic energy of bulk motion are much higher than the values from the numerical 
simulations even without heat conduction. 
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4.4. Direct observational implications 

Finally, we briefly discuss some direct observational implications of our models. Since 
for a correct construction of intensity maps one would need 3D models, the results of our 
2D calculations can only give a rough estimate of the observable intensities. Thus, we have 
calculated angle- averaged intensity profiles for Ha and soft X-ray emission based on the 
cylindrical symmetry in our models. 

The Ha emissivity is calculated according to the table (case B) in Osterbrock (1989). 
We use a T,7°-^^ fit for the temperature dependence and set the emissivity to zero where the 
degree of hydrogen ionization is below 10~^. We do not account for absorption. 

In Figure 22 we plot the angle-averaged Ha intensity profiles at two different times 
{t = 0.22 Myr, during which intense structure formation occurs in the H II region, and 
t — 3.365 Myr, at a late stage after the LBV phase of the star). For comparison, we have 
also plotted the Ha intensity profiles of the corresponding models without stellar wind for 
the same times. 

At t = 0.22 Myr the pure H ll region model without stellar wind shows the intensity 
profile of a spherical, emitting volume that has just started to expand. On the other hand, 
the intensity profile of the model with stellar wind shows a dominant peak between r — 8 
and 9 pc. This peak represents the global intensity maximum at that time and is produced 
by the emission of the dense shell swept up by the hot bubble of shocked stellar wind gas. A 
secondary, much smaller peak between r = 13 and 14 pc originates from the trapping of the 
ionization front in the dense outer shell fragments piled up by the expanding H ll region. 

Comparing these data with the intensity profiles a,tt — 3.365 Myr shows that the further 
expansion of the bubble and H ll region generally lowers the Ha surface brightness, as 
expected for comparable Lyman continuum fluxes. In particular, when the SWB completely 
overtakes the H II region and enlarges the whole structure, the Ha surface brightness is lower 
than in the pure H ll region model at the same time. 

Due to the facts that the photoionized region is more or less limited to the illuminated 
inner part of the shell and that the Ha emissivity in the hot cavity is very low, the intensity 
proflle of the combined SWB / H II region appears slightly limb brightened. The Ha emission 
from the WR bubble is barely visible above the Ha background from the MS bubble. It might 
be that in our simulation the shell of the WR bubble is too hot (e.g. due to the neglect of 
heat conduction) and thus emits more strongly in X-rays (see below) than in the optical. 

For the X-ray intensity proflles in the energy range from 0.5 keV to 3.0 keV we use 
the emissivity calculated with the Raymond & Smith (1977) program for cosmic chemi- 
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cal composition (Allen 1973). The emissivity is set to zero for temperatures below 10^ K; 
absorption is not considered. In Figure 23 we display the angle-averaged soft X-ray in- 
tensity profiles at t = 1.0 Myr, 3.30 Myr (before the LBV-phase), and 3.365 Myr (after the 
LBV-phase). The two snapshots before the LBV stage show intensities between 10^^'^ and 
10~^ergss~^ cm~^ sr~^ for lines of sight through the hot bubble. After the LBV stage, the 
soft X-ray intensity close to the center of the bubble is strongly enhanced by more than three 
orders of magnitude to values above 10~^ ergs s~^ cm~^ sr~^ when the shocked WR wind hits 
the LBV ejecta. The soft X-ray emission comes mainly from the inside of the shell of the 
swept-up LBV wind, before it breaks apart. 

The WR nebula RCW 58 around HD 96548 (= WR 40) is a possible candidate for an 
observed counterpart of the WR nebula that forms in our model calculation att — 3.365 Myr. 
The WR star is currently of type WN8 and there are indications that the star passed the 
LBV stage (GMLl; Humphreys 1991). The mean radius of the optical nebula is 3.5 pc (Chu 
1982) or 2.5 pc (Arthur et al. 1996) depending on the estimation of the distance to HD 
96548. Thus, the evolutionary state of the nebula should be approximately comparable to 
the nebula in our model. 

Though the strong increase of the surface brightness in our model after the LBV phase is 
basically in agreement with the fact that up to now only WR bubbles (and no MS bubbles) 
have been observed in X-rays, the comparison with RCW 58 shows that our model also 
suffers from the same problem that all analytical and numerical models thus far have: The 
X-ray luminosity is much higher than observed. In our model, the X-ray luminosity of the 
WR bubble in the energy range from 0.5 keV to 3.0 keV is 7 x 10^^ ergs s^^ whereas Moffat 
et al. (1982) report that no X-ray emission from RCW 58 was detected with the EINSTEIN 
satellite, meaning the X-ray luminosity of RCW 58 in this energy range must be below 
10^2 ergs s-^ (Arthur et al. 1996). 

4.5. The impact of the ambient medium, the stellar parameters, and 

unconsidered physics 

Due to the immense computational efforts we can only show results for one set of ambient 
medium parameters (no = 20 cm~^ and Tq — 200 K). We shall utilize the scaling laws for H ll 
regions in ionization/recombination and heating/cooling equilibrium with negligible gravity 
(see table 2 of Yorke et al. 1982) and briefly discuss the changes of the results that we expect 
if other values for the ambient medium are used. Separating the evolution of the H II region 
and SWB, we note that if the density in the ambient medium is modified by a factor 6 while 
assuming the same equilibrium temperature and same ionizing luminosity, all linear scales 
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A and time scales f transform according to A oc f oc velocities remain unchanged, 

and mass outflow rates and mechanical luminosity (as well as other rates of energy changes) 
transform as 5~^/^. Values for the various energies given at a particular time scale as 5~^. 
Thus, an increase in ambient density results in smaller H II regions, lower energy injection 
rates and proportionally lower values for the ionization, thermal, and kinetic energy. These 
same scaling laws can be derived in a less rigorous manner from our equations 2 to 5. 

Modifying the temperature of the ambient medium docs not strongly affect the evolution 
as long as the gas is not ionized and the pressure in the H II region is much higher than the 
pressure of the ambient medium. 

With respect to the SWB, equation 6 shows that the radius of the SWB is proportional 
to Po^^^, i.e. at comparable times it is smaller (larger) by 6^^^ in the case of an ambient 
density higher (lower) by a factor 5. In contrast to the H II region, the growth of the 
thermal and kinetic energy of the SWB does not depend on the ambient density, as long as 
the assumptions under that equations 9 and 10 were derived remain fulfilled (e.g. thermal 
pressure in the hot bubble is much higher than the thermal pressure in ambient medium). 

It is difficult to estimate the impact of a variation of the ambient conditions on the 
highly non-linear H II region / SWB interaction effects. Basically, an increase of the density 
in the ambient medium leads to both a smaller SWB and a thicker shell of swept up material. 
Thus, the ratio of shell thickness to SWB diameter grows which increases the unstable length 
scales and tends to inhibit thin shell instabifities with all their associated effects. On the 
other hand, dumpiness in the circumstellar medium may exist anyway, abolishing the need 
to trigger the interaction effects by shell instabilities. It is clear that numerical simulations 
under a variety of different ambient conditions are highly desirable in order to study these 
effects in more detail. They will become available with increasing computer power. 

It should also be noted here that the stellar parameters which wc use as time- dependent 
boundary conditions in our simulations are still somewhat uncertain. E.g. Martins, Schaerer, 
& Hillier (2002) showed that the determination of the effective temperature of O dwarfs based 
on non-LTE line blanketed atmosphere models including stellar winds results in differences 
to previous calibrations up to a few 1000 K. The mass-loss rate and its temporal variation 
in the short and violent evolutionary stages fike the LBV phase are probably even less well 
known. We have chosen the set of stellar parameters from GMLl because it is current state- 
of-the-art, sufficiently time-resolved, and we wanted to directly compare our 2D simulations 
to their ID calculations. 

Nevertheless, deviations of the stellar parameters from this set will certainly influence 
the results of the model calculations. Returning to the scahng laws discussed by Yorke et 
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al. (1982) and scaling the hydrogen-ionizing photon luminosity by I but keeping density 
and temperature constant, we find that lengths and times scale according to A oc f oc 
velocities remain unchanged, rates of energy changes scale like P^^, and values for energy 
injected into the ISM at a given time like /. In a similar manner, a higher stellar wind 
luminosity will increase the size of the SWB and the energy transferred to the shell at a 
given time. Note that although the analytical description of the SWB evolution depends on 
the stellar wind luminosity, i.e. only on the product of mass-loss rate and wind terminal 
velocity squared, this does not account for cooling in the hot bubble which will be enhanced 
when the mass-loss rate rises without a proportional increase in wind velocity. 

The impact of changes in the temporal evolution of the stellar parameters (e.g. a 
different number of outbursts of varying length) on the results of our calculations are difficult 
to estimate. Highly non-linear effects appear when e.g. a slow shell is overtaken by a fast 
wind and instabilities arise. Interaction processes between the SWB and the H ll region may 
also be amplified or damped by variations of the stellar parameters. Numerical simulations 
for different stellar evolutionary histories are necessary to study these effects in detail. 

Some physical processes are not yet considered in our models and we want to briefly 
discuss their possible influence here. One of these processes is heat conduction which may 
play a role for the energy transfer in the system. Though classical theory as well as our 
numerical models without heat conduction predict that the cooling time for the interior of 
the MS bubble is longer than the lifetime of the star, there is indirect evidence that, at 
least in NGC6888, the MS bubble is cold at the time when the WR wind blows out into 
the MS bubble (Mac Low 2000). A possible explanation is mass loading of the MS bubble: 
Heat conduction evaporates cold material from the shell or from immersed clouds which 
enhances the density in the bubble. The higher density reduces the coohng time of the 
bubble. Such evaporation from a cloud embedded in hot interstellar gas has recently been 
observed (Chu 2000). We will discuss the case of NGC6888 in a forthcoming paper that 
describes numerical simulations of the circumstellar medium around a 35 Mq star, a model 
which is more appropriate for NGC6888. 

On the other hand, magnetic flelds provide an efficient mechanism to reduce the effi- 
ciency of heat conduction. In a plasma the electrons are constrained to move along magnetic 
field lines. If the magnetic field lines are aligned with the shell or if they are tangled, the 
effective path lengths for fast electrons moving from the hot to the cold gas are much longer 
and heat conduction by electrons is less efficient. Saturation efi^ects of heat conduction due to 
electric fields arising from charge separation further complicates the situation. All these ef- 
fects are worthy to be studied in detail, but the simultaneous inclusion into multidimensional 
hydrodynamical simulations remains a computational challenge for the future. 
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We only allowed for isotropic winds in our models. However, line-driving theory for 

rotating stars suggests that mass outflow from the poles should be both more vigorous and 
faster. As already mentioned in section 2.4, Brighenti & D'Ercole (1997) and Frank et al. 
(1998) showed that anisotropic winds can indeed influence the formation and appearance of 
circumstellar shells, producing lobes or bipolar outflows. Similar effects may be expected for 
our models when anisotropic winds are employed, but details remain to be explored. 

Other physical effects not yet considered in our models include variations of the metal- 
licity in the circumstellar gas and the photodissociation of molecular hydrogen. In its later 
stages, the star is expected to eject gas that is enriched with metals. This enhancement 
of the metallicity in the circumstellar medium leads to stronger cooling and thus influences 
the circumstellar evolution. The dissociation of molecular hydrogen by stellar photons with 
energies below the Lyman threshold adds additional (exterior) layers to the basic structure 
shown in Figure 1. The H ll region evolves into the photodissociated region which itself 
expands into the molecular gas. Interaction processes between the photodissociated region 
and the H ll region similar to those between the H ll region and the SWB may arise. 

Finally, large scale density gradients such as those encountered at the edge of the clouds 
will give rise to a variety of additional dynamical effects, e.g. champagne flows of H ll regions 
and break-out of hot gas as well as the interactions between these two phenomena. 

5. Summary and conclusions 

Our numerical models for the evolution of circumstellar gas around a 60 star show 
that the interaction of the stellar wind bubble with the stellar radiation field can strongly 
influence the morphology of the circumstellar medium. The rearrangement of circumstellar 
gas by the stellar wind influences the way it reacts to the ionizing radiation. Density clumps 
formed in the shell of gas swept up by the stellar wind bubble cast shadows into the H II 
region. The resulting pressure gradients force material into the shadowed regions enhancing 
their density and forming neutral "spokes" , which subsequently raise the mass of the clumps 
in the shell. The H ll region is extended in directions free of clumps. Thus, the formation 
of these elongated H ll region "fingers" in our model occurs not only due to the ionization 
front instability described by Garcia-Segura & Franco (1996), but is triggered and amphfied 
by the redistribution of mass by the action of the stellar wind shell. 

These results also shed light on the open question whether the complex structures that 
can be found in H ll regions are primordial, i.e. relics from the time before the gas became 
ionized, or formed by dynamical processes in the course of the H ll region evolution. While 
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there are strong observational indications that the former plays an important role, our results 
give support to the idea that the latter cannot be completely excluded: Intense structure 
formation in H II regions with strong stellar winds can occur even if the neutral ambient 
medium was initially homogeneous. 

Nevertheless, these structures are only a temporary phenomenon because the extended 
H II region is eventually swept up by the stellar wind shell. If we compare our results at 
t = 3.30 Myr with the ID solution of GMLl, we see that, from a morphological point of view, 
there is basically little difference in the overall structure except for the appearance of an H II 
region at the inner part of the remnant shell and a few filaments in the hot bubble. Thus, 
the approach of GMLl to use the ID solution for the estimation of the initial conditions for 
closer studies of the LBV and subsequent WR stage appears to be vahd. 

The structure formation induced by the interaction of the stellar radiation field with the 
SWB also has implications for the energy balance of the circumstellar gas, mostly via the 
effect that denser gas has a shorter recombination time and a higher coohng efficiency. This 
can be seen in the decrease of ionization and thermal energy of warm photoionized gas in 
the circumstellar medium during the time when the structure formation occurs in the first 
few 10^ yr. Furthermore, the ionization energy lags behind the corresponding value in the 
model without wind by 0.1 — 0.4 dcx and the thermal energy of warm gas by 0.2 — 0.6 dcx 
for almost the entire lifetime of the star. This also implies that the H II regions around stars 
with winds can have a higher emission measure than undisturbed ones. 

The ionization energy and the thermal energy of warm gas in the circumstellar medium 
are lower compared to the models without wind. Though the total mechanical wind energy 
is negligible compared to the accumulated energy of the Lyman continuum photons, the 
simulation with a stellar wind results in a kinetic energy in the circumstellar medium which 
is 4 times higher than in the model without wind. The total thermal energy is almost 
twice as high (or, subtracting the initial thermal energy of the background gas, the thermal 
energy that is added to the system is enhanced by a factor of 2.4). The energy transfer 
efficiency of the stellar Lyman continuum radiation over long time scales is so low because 
this radiative energy is mostly used to maintain the photoionization of hydrogen; it is lost 
from the system when the hydrogen recombines into levels above the ground state. By 
contrast, most of the energy of the stellar wind is converted into thermal energy of hot gas 
that accelerates the shell. The kinetic energy can be accumulated in the system for a long 
time unless it is dissipated, and the thermal energy of hot gas can also be saved if the density 
in the bubble is sufficiently low. For a plasma of solar chemical composition with T = 10^ K 
and n — 10~^cm~^, the cooling time in coUisional ionization equilibrium is several 10^ yr. 
Although the LBV phase of the star induces very rapid changes in the energy balance of the 
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circumstellar medium, its impact over longer time scales is limited due to the brevity of the 
LBV eruption. 

The above conclusions regarding the inefficiency of converting the energy flux of pho- 
toionizing UV photons into kinetic energy are modified in the presence of large-scale density 
gradients (i.e. at the edge of a molecular cloud). Generally speaking, the resulting "cham- 
pagne fiow" can lead to an efficiency of about 1 % for converting the stellar UV flux into 
kinetic energy of expansion (see Yorke 1986). The exact value for efficiency depends on the 
details of the problem, however; it is conceivable in some cases that the champagne phase is 
negligibly short compared to the lifetime of the star. We note that simulations which assess 
the role of stellar winds in the presence of champagne flows are currently rare. Comeron 
(1997) examined the evolution of wind-driven H ll regions across a density discontinuity for 
various parameters but without the consideration of energy transfer efficiencies. 

For the simulations of HII region evolution without a stellar wind the analytical pre- 
diction for the kinetic energy in the system at the end of the star's lifetime differs from our 
numerical result by less than 30%. The analytical solution for the energy transfer efficiency 
in the case of the combined H ll region / SWB is just an upper limit and overestimates 
the transfer efficiencies into kinetic energy and thermal energy by factors of 6.7 and 9.7, 
respectively. 

This detailed examination is suitable to improve studies of the energization of the ISM 
in the solar neighborhood (e.g. Abbott 1982), which use energy transfer efficiencies as the- 
oretical input. It is also important when considering the heating of the galactic disk as a 
global phenomenon self-regulated by star formation. To address this the effects of overlap- 
ping H II regions and SWBs in OB clusters and associations still need to be assessed, the 
subject of a future investigation. 
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ambient interstellar gas 

Fig. 1. — Schematic structure of an SWB: Tgi marks the position of the reverse shock, Tc 
the contact discontinuity, rs2 the forward shock of the stellar wind bubble, ri the ionization 
front, and the forward shock of the H ll region expansion. 
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Fig. 2. — Timc-dcpcndent stellar parameters used as boundary conditions for the calculation 
of the 60 M0 model. Mass-loss rate (upper left diagram), terminal velocity of the wind (upper 
right), effective temperature (lower left), and luminosity (lower right). All parameters are 
adopted from GMLl. 



-39- 




0.0 0.2 0.4 0.6 0.8 1.0 

time [Myr] 

Fig. 3. — Resolution study for the 60 Mq model. is the total kinetic energy of bulk 
motion in the system, Et the thermal energy, and Ei the ionization energy (13.6 eV per 
ionized hydrogen atom). 
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Fig. 4. — Circumstellar mass density and velocity field (top panel), temperature (center 
panel), and degree of hydrogen ionization (lower panel) for the 60 model at age 10^ yr 
(high resolution run). The velocity arrows in the free- flowing wind zone have been omitted 
to prevent confusion. The star is located in the center of the coordinate system. Please note 
the different length scales. 




Fig. 5. — Same as Fig. 4 but at age 2 x 10^ yr. The velocity arrows in the free-flowing wind 
zone and in the hot bubble have been omitted to prevent confusion. 
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Fig. 7. — Same as Fig. 5 but at age 0.1 Myr. 
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Fig. 8 — Same as Fig. 5 but at age 0.22 Myr. 
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Fig. 10. — Same as Fig. 9 but at age 1 Myr. 




Fig. 11. — Same as Fig. 10 but for the medium resolution run. 
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Fig. 12. — Same as Fig. 11 but at age 3.30 Myr. Only mass density and velocity field (left 
panel) and temperature (right panel) are shown and the displayed area is enlarged once 
again. 
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Fig. 17. — The temporal evolution of kinetic, thermal, and ionization energy in the 60 Mq 
model. For details see text. 
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Fig. 18. — Ratio of energies in the 60 Mq model with and without a stellar wind. 
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Fig. 19. — Energy loss due to cooling in the 60 M© model. 
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Fig. 20. — Energy transfer efficiency witfi respect to tlie total energy input in the 60 
model without wind. 
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Fig. 21. — Energy transfer efficiency witfi respect to tlie total energy input in the 60 
model. 
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Fig. 22. — Angle-averaged Ha intensity for the 60 Mq model with and without a stellar wind 
are compared at two evolutionary times. 
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Fig. 23. — Angle- averaged soft (0.5 . . . 3.0keV) X-ray intensity for the 60 model with 
stellar wind at selected evolutionary times. 
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Tablc 1: Hydrodynamical models of SWBs / H ll regions around single early- type stars. 



Ref. 


Dimensions 


Geometry 


Linear size 


Resolution^ 








(pc) 


(10-3) 


1 


2 


cylindrical 


7.8... (3. 2 X 10-3) 


5 


2 


2 


spherical 


0.1 


5 


3 


1/2 


spherical 


70 ... 45 / 19 ... 3 


2.5 / 2.5 . . . 1.25 


4 


2 


spherical 


0.1... 2 


5 


5 


2 


spherical 


6.6 / 1 


1.3 


6 


2 


cylindrical 


0.2 


3.9 


7 


2 


cylindrical 


6.5 


2.5 


8 


2 


cylindrical 


64 


8.0... 0.125 


Ref. 


Duration 


Ionization 


Heating 


Cooling 




(Myr) 








1 


0.02... (1.5 X 10"^) 


thc^rmHl 


mechanical 


CIE"^ 


2 


10-3 


no 


mechanical 


isothermal EOS*^ 


3 


4...5 / (28... 7.2) X 10-3 


no 


mechanical 


CIE + cutoff 


4 


0.02... 0.33 


PIE'^ 


mcch. + PIE 


CIE + cutoff 


5 


0.023 


no 


mechanical 


CIE 


6 


8 X 10-^ 


no 


mechanical 


CIE 


7 


0.015 


no 


mechanical 


CIE + cutoff 


8 


4 


time dependent 


mech. + radiative 


exphcit or CIE* 


Ref. 


Wind asymmetry 


iw(t) 


J^Lyc(i) 


no{r) 


1 


no 


fixed / two-level 





const. / composite 


2 


no 


two-level 





const. 


3 


no 


variable 


variable 


const. / from ID 


4 


no 





fixed 


const. / power law 


5 


yes 


two-level with trans. 





pre WR 


6 


yes 


three-level 





const. 


7 


no 


const. 





const. 


8 


no 


variable 


variable 


const. 



"in units of linear size 

''cooling function based on the assumption of collisional ionization equilibrium 
■^equation of state 

''photo ionization equilibrium calculated assuming that the gas is fully ionized inside the H n region 

"^thermal and radiative ionization of hydrogen 

•^explicit calculation of important cooling processes for lower temperatures and CIE for high temperatures 

References. — (1) Rozyczka (1985), Rozyczka & Tenorio-Taglc (1985a,b); (2) Stone et al. (1995); (3) 
GMLl, Garci'a-Segura et al. (1996b); (4) Garci'a-Segura & Franco (1996); (5) Brighenti & D'Ercole (1997); 
(6) Frank et al. (1998); (7) Strickland & Stevens (1998); (8) This work. 
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Table 2: The energy components at the end of the simulation. 



Model 






-E't.cold 


-^t,Wcirm 


-E't.hot 




(lO^^ej-gs) 


(lO^^ej-gs) 


(lO^^ergs) 


(lO^^ergs) 


(lO^^ergs) 


windless 


3.5 


13 


5.5 


4.7 





with wind 


14 


10 


7.9 


2.1 


9.7 



Note. — The thermal energy of the cold component, -Et,coid) contains the internal energy of the initially 
unperturbed ambient medium (3.4 x 10^^ ergs) that has to be subtracted whenever the input of thermal 
energy into the system is considered. 
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Table 3: The energy transfer efficiencies at the end of the simulation. 



Model 


£k 








(10-4) 


(10-4) 


(10-4) 


windless 


3.3 


12 


6.4 


with wind 


13 


9.1 


15 



